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This  work  presents  a  circuit  model  for  calculating  the  total  energy  dissipated  into  neu¬ 
tral  species  for  nanosecond  pulsed  direct  current  (DC)  dielectric  barrier  discharge  (DBD) 
plasmas.  Based  on  experimental  observations,  it  is  assumed  that  the  nanosecond  pulsed 
DBD’s  which  have  been  proposed  for  aerodynamic  flow  control  can  be  approximated  by 
two  independent  regions  of  homogeneous  electric  field.  An  equivalent  circuit  model  is  de¬ 
veloped  for  both  homogeneous  regions  based  on  a  combination  of  a  resistor,  capacitors,  and 
a  zener  diode.  Instead  of  fitting  the  resistance  to  an  experimental  data  set,  a  formula  is 
established  for  approximating  the  resistance  by  modeling  plasmas  as  a  conductor  with  DC 
voltage  applied  to  it.  Various  assumptions  are  then  applied  to  the  governing  Boltzmann 
equation  to  approximate  electrical  conductivity  values  for  weakly  ionized  plasmas.  The 
developed  model  is  then  validated  with  experimental  data  of  the  total  power  dissipated  by 
plasmas. 


Nomenclature 


B  Magnetic  Field,  T 

C  Capacitance,  F 

E  Electric  Field,  V/m 

J  Current  Density,  A/m2 

l  Chordwise  Length,  m 

m  Mass,  kg 

Na  Number  Density  of  Atoms,  m-3 
Ne  Number  Density  of  Electrons,  m-3 
P  Power,  W 

Q  Energy,  J 

R  Resistance,  B 

v  Velocity,  m/s 

V  Voltage,  V 

Vapp  Applied  Voltage,  V 

Vvoi  Volume,  m3 

we  Drift  Velocity  of  Electrons,  m/s 

Te  Electron  Flux,  m-2s-1 

ea  Relative  Dielectric  Constant  of  Air 
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I.  Introduction 

The  need  for  improved  control  over  aerodynamic  flow  separation  has  increased  interest  in  the  potential  use 
of  plasma  actuators.  The  inherent  advantages  of  plasma  actuator  flow  control  devices  include:  fast  response 
time,  surface  compliance,  lack  of  moving  parts,  and  inexpensiveness.  However,  it  has  been  established  that 
the  actuators  which  affect  the  flow  via  directed  momentum  transfer  are  not  effective  at  Mach  numbers  as¬ 
sociated  with  most  subsonic  aircraft  applications.  Recently,  Roupassov  et.  al10  demonstrated  that  pulsed 
plasma  actuators,  in  which  energy  imparted  to  the  flow  appears  to  effectively  control  flow  separation,  seem 
to  be  suitable  at  Mach  numbers  (M«0.3)  beyond  the  capabilities  of  the  current  plasma  induced  momentum 
based  approaches. 

Given  the  fundamental  differences  between  the  novel  pulsed  discharge  approach  and  the  more  conventional 
momentum  based  approaches,  there  is  a  need  to  develop  an  effective  and  efficient  model  for  the  energy 
delivered  to  the  flow  by  the  plasma.  Once  calculated,  that  value  can  be  input  to  a  computational  fluid 
dynamics  solver  as  an  energy  source  term  resulting  in  a  coupled  fluid/plasma  dynamics  model.  Multiphysics 
models  of  this  type  are  required  in  order  to  study  detailed  flow  characteristics.  However,  detailed  numerical 
simulations  involving  plasma  kinetics  are  computationally  prohibitive  for  a  variety  of  coupled  fluid/plasma 
design  problems.  To  address  this  issue,  efficient  circuit  element  models  have  been  introduced  to  approximate 
the  complex  processes  within  plasmas.  Circuit  models  such  as  those  by  Orlov  et.  al6  rely  on  empirical 
constants  which  may  not  be  applicable  to  nanosecond  pulsed  discharges.  To  date,  an  approximate  model  of 
nanosecond  pulsed  plasma  actuators  has  not  been  developed.  This  paper  deals  primarily  with  establishing  a 
flexible  model  with  relevant  physics  that  could  be  implemented  as  an  approximation  for  the  energy  dissipated 
within  a  plasma  for  any  pulsed  DC  DBD  configuration.  Among  the  other  goals  in  this  paper  is  to  probe  into 
the  background  processes  that  occur  within  plasmas  and  incorporate  that  knowledge  into  the  model. 

II.  Lumped  Element  Circuit  Model 

One  of  the  primary  assumptions  in  creating  this  model  is  that  nanosecond  pulsed  DBD’s  can  be  approximated 
by  two  independent  regions  of  homogeneous  electric  field.  One  such  region,  dubbed  the  ‘hot  spot’  is  the 
region  adjacent  to  the  powered  electrode.  This  region  makes  up  a  small  portion  of  the  total  discharge  area 
but  was  observed  to  be  an  important  component  of  the  plasma  discharge  and  necessary  to  obtain  agreement 
with  experimentally  measured  shock  wave  dynamics  by  Roupassov  et.  al.10  The  other  region,  dubbed  the 
‘tail,’  encompasses  the  rest  of  the  plasma  discharge  and  extends  to  the  edge  of  the  dielectric.  As  both  regions 
are  independent,  the  model  presented  in  this  paper  consists  of  a  single  network  for  each  region  containing  a 
resistor,  capacitors,  and  a  diode. 

As  shown  in  Fig.  1,  circuit  elements  that  were  used  to  model  the  plasma  include:  an  air  capacitor  Ca,  a 
dielectric  capacitor  Cd,  a  resistor  Rf,  and  a  zener  diode  Df.  The  air  capacitor  represents  the  capacitance 
between  the  dielectric  surface  and  the  exposed  electrode.  The  dielectric  capacitor  represents  the  capacitance 
between  the  dielectric  surface  and  insulated  electrode  and  is  proportional  to  the  thickness  of  the  dielectric 
layer.  Thus  the  dielectric  layer  in  the  form  of  both  its  thickness  and  the  value  of  its  dielectric  constant 
plays  an  important  role  in  determining  the  effectiveness  of  the  plasma  actuator.  Finally,  the  zener  diode, 
introduced  by  Orlov  et.  al,6  is  utilized  in  the  model  to  enforce  an  energy  threshold  value  below  which  plasma 
will  not  form. 

Since  a  uniform  charge  distribution  along  the  top  of  the  dielectric  is  assumed,  the  typical  asymmetric  2-D 
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Figure  1.  Electric  circuit  model  of  a  dielectric  aerodynamic  plasma  actuator. 


plasma  actuator  geometry  featured  in  Fig.  1  can  be  simplified  to  a  series  of  homogeneous  symmetric  regions. 
As  shown  in  Fig.  2,  these  regions  include:  an  anode  sheath,  ’hot  spot’,  and  ‘tail.’ 


Exposed  Electrode 


"Virtual"  Electrode 


Figure  2.  Plasma  discharge  regions. 


This  assumption  results  in  a  series  of  coupled  1-D  models  that  account  for  the  chordwise  variation  along  the 
actuator.  Fig.  3  shows  the  simplified  circuit  model  within  each  homogeneous  region. 


A.  Circuit 

As  displayed  in  Fig.  1,  the  lumped  element  circuit  is  a  function  of  the  two  capacitance  values,  Ca  and  Cd- 
In  this  model,  the  air  is  treated  as  both  a  conductor  to  generate  a  physical  relationship  for  the  resistance 
Rf  and  a  parallel  plate  capacitor  to  generate  Ca -  An  advantage  of  modeling  the  plasma  as  a  conductor  in 
addition  to  a  parallel  plate  capacitor  is  that  it  generates  a  physical  relationship  for  the  resistance,  Rf ,  a 
value  that  is  traditionally  empirically  determined.  The  air  gap  capacitor  can  be  modeled  as4 


Ca  = 


K 


(i) 


where  Aa  is  the  cross-sectional  area  of  the  air  and  ha  is  the  approximate  height  of  the  plasma  region  of 
interest.  The  height  of  the  plasma  has  been  shown  by  Roupassov  et.  al10  to  be  approximately  independent 
of  applied  voltage  for  nanosecond  pulsed  DBD  actuators.  As  displayed  in  Fig.  4,  Aa  is  the  product  of  the 
spanwise  length  of  the  actuator  za,  and  la  is  the  chordwise  distance  from  the  exposed  electrode  to  the  end 
of  the  dielectric  region. 


The  capacitive  element  corresponding  to  the  dielectric  can  be  modeled  as4 

eo  edAd 


Cd  = 


hd 


(2) 
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Figure  3.  Region  of  homogeneous  potential. 


Figure  4.  Sketch  of  the  capacitive  air  element. 


where  A d  is  the  cross-sectional  area  of  the  dielectric  capacitive  element  and  h d  is  the  height  of  the  dielectric 
barrier  layer.  As  displayed  in  Fig.  5,  Ad  is  the  product  of  the  spanwise  length  of  the  actuator  za  and 
the  width  of  the  dielectric  region.  Treating  the  plasma  as  a  conductor,  the  resistance  for  DC  voltage  is 
proportional  to  crp,  Aa,  and  ha  and  can  be  given  as4 


Rf  = 


ha 


(3) 


Starting  from  Kirchoff’s  circuit  laws,4  the  governing  differential  equation  for  the  voltage  drop  experienced 
by  the  air  gap,  Ay,  is  given  by 


dAV(t)  =  dVapp  f  Ca  \  A V(t) 

dt  dt  \Ca  +  Cd  )  KRs(t){Ca  +  Cdy 

^  _  f  1  if  |E|  >  Ecrit , 
lb  if  |E|  <  Ecrit. 


(4) 

(5) 


where  Vapp  is  the  applied  voltage  and 
magnitude,  given  as 


k  is  the  contribution  from  the  zener  diode. 


If  the  electric  field 


(6) 


is  greater  than  the  breakdown  electric  field,6  Ecr^,  required  to  ionize  air,  then  k  takes  on  a  value  of  one, 
otherwise  it  is  zero  to  signify  that  plasma  has  not  formed. 
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Figure  5.  Sketch  of  the  dielectric  capacitive  element. 


B.  Conductivity 

To  effectively  calculate  the  resistance  governed  by  Eq.  3,  an  expression  must  first  be  developed  for  the 
electrical  conductivity,  crp,  of  the  plasma.  This  value  is  one  that  traditionally  requires  a  numerical  approach. 
To  simplify  the  problem  to  a  point  where  an  analytic  formulation  can  be  used,  numerous  simplifying  as¬ 
sumptions  were  used  and  are  described  in  the  following  paragraphs. 

For  any  plasma,  the  resulting  electric  current  is  composed  of  two  primary  terms:  the  current  from  electrons 
and  that  from  ions.  As  the  drift  velocity  of  electrons,  we,  in  a  non-equilibrium  plasma  is  significantly  higher 
than  ions,  the  current  density  can  be  approximated  as  the  portion  from  electrons  if  the  number  densities, 
Ne  and  A^,  are  approximately  the  same.  Using  a  form  of  the  generalized  Ohm’s  Law,  the  current  density 
vector,  J  and  ap  respectively  can  be  written  as 

J  «  —eNe we  =  crpE,  (7) 

(jp  =  e  ( Nefjje  +  Nif^i) ,  (8) 

where  fa  and  fie  represent  the  ion  and  electron  mobilities  respectively.  Much  like  Eq.  7,  the  electrical  conduc¬ 
tivity  relation  can  be  simplified  using  the  concept  of  quasineutrality  which  is  defined  as  having  approximate 
equal  number  densities  for  charged  particles  of  opposite  polarity.  Thus  as  fie  is  typically  three  orders  of 
magnitude  larger  than  fa ,  it  is  a  good  assumption  to  approximate  the  electrical  conductivity  as  only  coming 
from  electrons  as  long  as  Ne  is  at  least  of  the  same  order  of  magnitude  as  A^.7  Quasineutrality  is  a  typical 
assumption  that  is  valid  as  long  as  the  plasma  being  modeling  is  far  enough  away  from  the  powered  electrode 
to  avoid  the  boundary  layer  in  plasma  physics  called  the  sheath. 

Since  a  pulsed  DC  voltage  is  assumed,  the  activation  of  the  external  electric  field  will  follow  the  voltage 
waveform  as  a  step  function.  Thus  two  expressions  will  be  required  for  crp,  where  the  first  is  valid  for 
the  period  when  an  external  electric  field  is  applied,  as  shown  in  Fig.  9  from  0-65  ns,  and  the  second 
when  the  voltage  drop  over  the  air  gap  is  zero.  For  the  portions  of  the  voltage  waveform  that  AU  is  zero, 
the  power  is  also  zero  according  to  Ohm’s  Law  and  thus  the  conductivity  during  this  time  is  of  no  importance. 

To  generate  a  analytic  formula  for  the  electrical  conductivity,  a  distribution  function  must  be  introduced  to 
describe  the  physical  evolution  in  the  number  of  particles,  /(v,  J,  r,  £),  defined  such  that  /(v,  J,  r,  t )  dv  is  the 
number  of  particles  in  a  unit  volume  located  at  point  r,  time  £,  internal  quantum  number  J,  and  differential 
velocity  range  v  +  dv.  Using  this  distribution  function,  the  number  of  particles  at  point  r  and  time  t  can  be 
defined  as 

N(r’t)  =  '52  [ /(vVM)dv.  (9) 

This  distribution  function  allows  a  mathematical  description  to  be  developed  for  the  temporal  evolution  in 
the  number  of  particles  resulting  from  particle  collisions  within  a  control  volume.  The  time  rate  of  change 
in  the  number  of  particles  due  to  externally  applied  fields  can  be  described  as11 

Df  =  /(v  +  dv,  J,  r  +  dr,  t  +  dt)  -  /(v,  J,  r,  t) 

Dt  dt 
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A  partial  differential  equation  can  be  developed  to  describe  Eq.  10  using  formulas  for  the  time  rate  of  change 
of  v  and  r  established  by  using  the  equation  of  motion  in  the  form  of 


F 

m  ’ 


dv 
dt 

dr 


and  the  chain  rule  of  calculus  to  obtain  the  Boltzmann  kinetic  energy  equation  given  as 

Rl  =  ^l  ,  m9f  F  a/ 

Dt  dt  V  dr  m  dv' 


(11) 

(12) 

(13) 


In  order  to  approximate  the  total  derivative,  a  relaxation  time  can  be  introduced  defined  as  the  time  taken 
for  the  system  to  be  reduced  to  an  equilibrium  distribution  function.  The  tau  approximation  can  be  given 
as  r  ^  (Na<jea |v|)-1  where  aea  is  the  collision  cross  section  between  electrons  and  atoms,  |v|  is  the  average 
collision  velocity,  and  Na  is  the  number  density  of  atoms.11  If  the  number  of  particles  within  a  control 
volume  is  defined  as  a  equilibrium  distribution  function,  /o,  when  each  particle  is  at  the  same  energy  level 
as  its  neighbor,  then  the  particle  evolution  over  time  due  to  pairwise  collisions  after  an  external  force  has 
been  applied  can  be  given  as11 

^  =  -  — ■  (14) 

Dt  T 

As  r  only  accounts  for  collisions  between  electrons  and  neutral  atoms,  it  is  only  accurate  in  the  event  of 
a  weakly  ionized  plasma.  Plasma  actuators  considered  in  this  paper  traditionally  feature  a  low  degree  of 
ionization,  or  simply  the  amount  of  air  that  is  ionized,  and  thus  can  be  treated  in  a  weakly  ionized  limit. 
In  terms  of  the  momentum-transfer  collision  frequency  which  can  be  defined  as  the  mass  corrected  rate  at 
which  a  particle  of  a  specific  species  collides  with  another,  the  criteria  for  a  weakly  ionized  plasma  can  be 
given  as7 

Vei  ^en*  (15) 

This  equation  requires  that  the  collision  frequency  between  electrons  and  ions  be  much  less  that  those  be¬ 
tween  electrons  and  neutrals.  Thus,  if  this  requirement  is  met,  the  collisional  occurrences  between  electrons 
and  charged  particles  can  be  effectively  ignored  and  the  relaxation  time  established  is  a  good  approximation 
of  the  total  time  rate  of  change  in  the  number  of  particles  within  a  control  volume. 


After  introducing  a  relaxation  time  to  approximate  Eq.  13,  an  equation  of  motion  describing  the  average 
velocity  of  electrons  can  be  established  by  multiplying  by  mev  and  integrating  over  the  electron  velocity.  An 
analytic  formulation  for  the  average  velocity  an  electron  experiences  due  to  an  externally  applied  field,  we, 
can  be  obtained  by  assuming  that  the  force  term  can  be  approximated  as  the  Lorentz  force, 

F  =  -eE  -  -  (v  x  B) ,  (16) 

c 

where  B  represents  the  magnetic  field  and  E  represents  the  electric  field.  As  plasma  actuators  have  no  applied 
magnetic  field,  B  can  be  set  equal  to  zero.  Therefore  the  equation  of  motion  for  an  electron  describing  we 
can  be  given  as 


dwe(t)  W e(t)  _  . 

me — -=—^-  +me - =  — eE  (t). 

dt  t 

Solving  Eq.  17,  a  linear  first-order  differential  equation,  an  integral  equation  is  obtained: 

[  E(t )  exp  (  —  )  dt. 


we(t *)  = - exp 

mP 


(17) 


(18) 


Eq.  18  can  be  solved  in  conjunction  with  Eq.  7  to  obtain  an  expression  for  the  time  varying  conductivity 
and  in  conjunction  with  Eq.  8  to  obtain  an  approximation  for  the  time  varying  electron  mobility  if  the  ion 
mobility  is  neglected.  The  electrical  conductivity  and  electron  mobility  respectively  can  be  given  as 


CTp(f)  = 


Ne{t*)e 


meE(t 

f 

meE(t*)  J0 


exp 


E(t)  exp 


t-t* 
N*cr*  v* 

a  u  eau 

t~t * 
N*<j *  v* 

au  eaV 


dt, 


dt , 


(19) 

(20) 
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where  aea  is  a  function  of  electron  energy  and  can  be  obtained  for  various  molecules  found  in  air  from 
Phelps.8  Numerical  values  used  in  the  model  are  included  in  the  Appendix.  The  electron  velocity  can  be 
obtained  by  assuming  a  Maxwellian  velocity  profile.  As  a  collection  of  electrons  within  plasma  have  a  range 
of  velocities,  the  Maxwellian  velocity  profile  represents  the  most  probable  distribution  of  these  velocities. 
Thus  the  distribution  of  velocities,  f(u,v,w),  can  be  given  by7 


f  (u,v,w)  =  As  exp 


\ \m  (u2  +  v2  +  w2) 
kbTe 


As  =  NP 


m 


27 rkbTe 


(21) 

(22) 


Using  the  non-relativist ic  definition  of  kinetic  energy,  a  relationship  can  be  established  for  the  kinetic  energy 
of  an  electron  that  is  valid  in  the  limit  |v|  <<  c,  where  c  is  the  speed  of  light.  Using  this  approximation,  the 
relativistic  formulation  of  the  kinetic  energy  can  be  approximated  as 

E„  =  mc2  -  me*  ee  imv*.  (23) 

V'l  -  PIT?  2 

Averaging  Eq.  21  and  using  the  non-relativistic  definition  of  kinetic  energy,  the  mean  kinetic  energy  of 
electrons  becomes7 

As\me  ( u 2  +  v2  +  w2)  exp  [—  \me  ( u 2  +  v2  +  w2)  / kbTe  1  du  dv  dw  3 

fffooV  , - r  X  ,  T  A  23  ,  '  ,  -  =  M,  (24) 

JJJ_oo  As  exp  [—  |me  [uz  +  vz  +  wz)  /  kbTe\  dudvdw  2 

where  kb  is  Boltzmann’s  constant.  From  the  definition  of  kinetic  energy,  the  relationship  between  Eav  and  |v| 
with  vector  components  (u,v,w)  can  be  established  and  the  average  thermal  velocity  of  electrons  becomes 


v  = 


/ 


3  kbTe 
me 


(25) 


The  required  inputs  for  Eqs.  19-20  include:  Ne  the  number  density  of  electrons,  Na  the  number  density  of 
atoms,  and  Te  the  temperature  of  the  electrons.  Among  these  values,  Na  can  be  assumed  to  be  constant  in 
time  as  the  number  density  of  atoms  is  significantly  higher  than  that  of  free  electrons. 


Many  other  models  incorporate  a  constant  electron  temperature  into  their  model.5,6  Using  experimentally 
measured  values  of  reduced  electric  field  strength,1  E/Na  vs.  electron  temperature  as  detailed  in  Fig.  6,  this 
model  calculates  a  new  electron  temperature  at  each  time  step  by  comparing  the  E/N  value  experienced  by 
the  plasma,  produced  from  using  Eqs.  4-6,  with  Fig.  6.  It  is  assumed  in  this  paper  that  the  effects  of  an 
applied  electric  field  have  an  instantaneous,  or  on  a  time  scale  much  faster  than  10-9  s,  effect  on  electrons. 


A  time- varying  differential  equation  that  governs  Ne  can  be  obtained  from  the  drift-diffusion  equations.  The 
electron  continuity  equation  that  governs  Ne  can  be  given  by 

^+V-re  =  a|re|-/3nine,  (26) 

where  Te  is  called  the  charged  species  flux,  a  is  the  Townsend  coefficient  of  ionization,  and  /?  is  the  re¬ 
combination  coefficient  between  electron  and  neutral  atoms.  By  simplifying  the  plasma  discharge  into  a 
combination  of  two  homogeneous  regions  plus  an  anode  sheath  and  by  invoking  the  irrotational  property  of 
electric  fields,  Eq.  26  can  be  simplified  by  ignoring  any  spatial  variation  in  the  number  density  of  electrons 
i.e.  re  becomes 

re  «  Nefie\E\.  (27) 

Assuming  the  number  densities  of  electrons  and  ions  are  equal  i.e.  in  the  quasineutral  region,  the  electron 
continuity  equation  for  air,  with  a  composition  of  80%  N2  and  20%  O2,  can  be  written  as 


dNe 

dt 


cw|iVeMeE|  -  0.80/3jV2-ZVg  -  O.2OP02NI 


(28) 
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14 


Figure  6.  Plot  of  electron  temperature  vs.  reduced  electric  field. 


where  aair  and  /?  respectively  can  be  approximated  as  having  the  form9 


(An 


Apex p  ( 


PN2  =  2.8  x  10  7  [cm3/s], 

P02  =  2  X  1CT7  ^^^)  [cm3/s], 


(29) 

(30) 

(31) 


where  A  and  B  are  empirical  constants  that  have  tabulated  values  of  15  and  365  respectively  for  air  at 
atmospheric  pressure.9 


C.  Discharge  Development 

The  final  remaining  unknown  required  to  close  the  equation  system  is  to  determine  how  the  electric  potential 
changes  over  the  horizontal  length  of  the  actuator.  To  accurately  model  this,  it  is  important  to  incorporate 
the  wall  effects  of  the  plasma  actuator.  In  terms  of  potential  variation,  these  wall  effects  attract  charged 
particles  of  opposite  polarity  and  shield  charged  particles  of  the  same  polarity.  Therefore  for  regions  beside 
the  anode  (powered  electrode  for  positive  pulses),  an  anode  sheath  is  developed  where  an  attraction  of 
electrons  occurs  and  a  repulsion  of  positive  ions  occur.  For  regions  above  the  dielectric  region  on  top  of 
the  grounded  electrode,  a  cathode  sheath  is  developed  where  positive  ions  are  collected  and  electrons  are 
repelled.  Fig.  7  illustrates  the  two  predominate  sheaths  that  are  developed  in  asymmetric  plasma  actuators. 

It  is  important  to  consider  the  effects  of  charge  collection  and  repulsion  in  these  regions  as  such  phenomenon 
can  have  large  effects  on  the  variation  of  the  electric  potential  over  the  chordwise  length  and  height  of  the 
plasma  discharge. 

Wall  Effects 

As  this  model  is  interested  in  solving  the  amount  of  energy  the  plasma  transfers  to  neutral  species,  the 
relative  importance  of  both  the  anode  and  cathode  sheath  regions  needs  to  be  established.  Energy,  Q,  is  a 
function  of  electrical  conductivity  and  electric  field  strength  and  can  be  given  as7 

Q  =  ap  \E2\Vvol.  (32) 
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Figure  7.  Collective  wall  effects  of  the  exposed  powered  electrode  and  virtual  electrode  (dielectric  region).  1- 
Anode  Sheath,  2-  Cathode  Sheath.3 


Eq.  32  shows  that  the  energy  transfer  from  both  sheath  regions  is  proportional  to  the  conductivity  of  the 
plasma  within  each  region.  As  the  cathode  sheath  region  typically  has  very  few  electrons  due  to  repulsion 
effects,  the  current  within  this  region  will  be  carried  by  positively  charged  ions.  The  mobility  of  such  ions 
are  significantly  less  than  that  of  an  electron,  so  as  a  first  order  approximation  if  the  local  charged  species 
number  density  and  electric  field  strength  are  assumed  equal,  the  conductivity  in  the  cathode  sheath  will  be 
significantly  less  than  the  conductivity  in  a  bulk  plasma  and  the  anode  sheath  i.e. 

(J S  —  c  ~  &S  —  CL'  (33) 

When  combined  with  the  fact  that  the  volume  of  both  sheath  regions  are  orders  of  magnitude  less  than  the 
bulk  plasma,  an  order  of  magnitude  approximation  to  the  plasma  discharge  can  be  obtained  by  ignoring 
the  cathode  sheath’s  effects  for  nanosecond  pulsed  plasma  discharges.  Although  the  cathode  sheath  is  ap¬ 
proximately  negligible  in  terms  of  energy  transfer,  the  higher  electric  field  strength  and  higher  conductivity 
present  in  the  anode  sheath  contains  important  physics  necessary  for  capturing  the  energy  transferred  by  a 
plasma  discharge. 

Fig.  2  shows  the  updated  ‘hot  spot’  and  ‘tail’  regions  with  the  anode  sheath  included  in  the  hot  spot  region 
adjacent  to  the  powered  electrode. 


Electric  Potential  Variation 


To  establish  a  variation  in  the  electric  potential,  the  governing  Maxwell  equations  can  be  used  which  are 
written  as4 


V  •  E  =  — , 

eo 

V  •  B  =  0, 

„  „  dB 

VXE=-^’ 

„  „  dE 

V  x  B  —  /i0J  + 


(34) 

(35) 

(36) 

(37) 


in  differential  form  where  pf  =  e  (rii  —  ne)  is  the  net  charge  density  and  po  is  the  permeability  of  free  space. 
In  the  absence  of  a  time- varying  magnetic  field,  Eq.  36,  Faraday’s  Law  of  Induction,  simplifies  to 


V  x  E  =  0, 


(38) 


and  since  the  curl  of  E  is  zero,  the  electric  field  can  be  solved  for  as  a  potential  function  (j)  and  substituted 
into  Gauss’  Law,  Eq.  34.  The  resulting  equations  respectively  can  be  given  by 


E  =  -V0,  (39) 

V2<^=“-  (40) 

The  net  charge  density  within  the  quasineutral  region  of  a  plasma  is  equal  to  zero  as  ne  =  ni.  For  the  anode 
sheath,  the  net  charge  density  can  be  approximated  if  the  plasma  is  assumed  to  uniformly  distribute  its 
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charge  density.  Therefore  the  quasineutral  region  of  the  ‘hot  spot’  and  ‘tail’  feature  equal  number  densities 
of  both  electrons  and  ions  (ne  =  rii)  while  the  cathode  and  anode  sheaths  only  have  electron  and  ion  densities 
respectively.  If  this  is  assumed  then  the  anode  sheath  can  be  approximated  as  having  a  net  charge  density 
equal  to  the  quasineutral  region’s  calculated  electron  number  density.  By  assuming  the  anode  sheath  is 
devoid  of  any  positive  ions  and  has  a  time-dependent  number  density  of  electrons,  Poisson’s  equation  for  the 
sheath  can  be  given  by 

<P<t>  eNe(t) 

dx 2  e0  '  [  ’ 

If  this  form  is  assumed  for  a  single  Debye  length  (A d),  then  the  remaining  discharge  (rest  of  ‘hot  spot’  and 
‘tail’)  region  is  part  of  the  quasineutral  bulk  plasma  and  can  be  given  by  Laplace’s  equation 


(42) 


The  ordinary  differential  equations  for  electric  potential  variation  in  the  hot  spot  and  tail  regions  respectively 
can  be  solved  to  provide  an  approximate  1-D  spatial  variation. 


4>{x)  = 


^x2  +  ClX  +  C2 

C$x  +  C4 


if  x  <  \d 
if  x  >  Xd 


(43) 


Eq.  43  requires  a  total  of  4  boundary  conditions.  Those  can  be  summarized  as 


0  l(x  —  0)  —  kapp, 

(44) 

02  (^  =  A)  =  Vbreak 5 

(45) 

ll 

V 

II 

to 

II 

V 

(46) 

II 

<M  I  , 

^  1^ 
11 

11 

1^ 

(47) 

where  0i  is  the  potential  in  the  anode  sheath,  02  is  the  potential  in  the  quasineutral  region,  L  is  length  of 
the  actuator  and  Xd  is  a  Debye  length.  The  first  boundary  condition  is  Vapp  at  the  cathode  and  the  second 
is  based  on  experimental  observations  by  Roupassov  et.  al.10  It  was  observed  that  a  plasma  discharge  could 
be  approximated  as  stopping  at  the  edge  of  the  grounded  electrode  for  asymmetric  actuators  independent  of 
the  applied  voltage;  so  the  edge  represents  the  absolute  limit  of  ionization  or  the  breakdown  voltage  of  air. 
Using  Eqs.  44-47  as  boundary  conditions,  the  potential  becomes  Eq.  48. 


eNe  2 
2e0 


eNe 

2cqL 


x0  -  + 


E>reafc  Ex  5 


x  +  V0 


app 


eNe_  2  _  eN_e  E  reak  Eipp 

2e0 LX0  eo  ^0  w  L 


(x  L )  +  Vbreak 


if  x  <  Xd 
if  x  >  Xd 


(48) 


The  Debye  length  provides  an  order  of  magnitude  approximation  for  the  extent  of  a  plasma  sheath  by  assum¬ 
ing  an  exponential  Boltzmann  distribution  in  the  charge  density  within  the  plasma  discharge.  Substituting 
this  into  Poisson’s  Eq., 


d 20  M 

e°^= eN° 


exp 


e0 

hTe 


-  1 


(49) 


where  N ^  is  the  charged  particle  density  far  away  from  the  electrode.  Taking  a  first-order  Taylor  expansion, 
the  Debye  length  can  be  given  as7 


A  d 


f  e0kbTe  \  2 
WooeV 


(50) 


Although  at  high  voltages,  a  first-order  approximation  fails,  Eq.  50  still  provides  an  order  of  magnitude 
approximation  of  the  extent  of  the  anode  sheath. 
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D.  Numerical  Procedure 

When  solving  for  the  energy  imparted  to  neutral  species,  a  coupled  equation  system  results  from  Eqs.  4  and 
28.  The  coupled  terms  include: 


Rf  oc  Ne, 

(51) 

re  oc  AV, 

(52) 

/?  oc  AV, 

(53) 

4>1  OC  Ne. 

(54) 

To  solve  the  resulting  equation  system,  the  Dormand-Prince  Runge-Kutta  method  was  employed.  This 
method  provides  an  efficient  way  to  incorporate  an  adaptive  step  size  that  is  important  for  computational 
efficiency  in  a  problem  that  requires  small  time  steps  for  convergence.  The  benefit  of  such  a  procedure  can 
be  illustrated  through  a  simplistic  example.  If  the  error  of  each  time  step  is  defined  as 

4  =  y(t)  -  yh 

(55) 

then  if  two  step  sizes  are  considered,  hi  and  /12,  the  error  of  each  iteration  and  their  relative 
can  be  given  as 

error  respectively 

V(t)  ~  Vnl  =  enl  =  ah 1» 

(56) 

y(t)  -  Vnl  =  en2  =  ah2, 

(57) 

3  ^ — N 
to  13 

1 

II 

© 

1 

to 

(58) 

Therefore  for  a  given  error  tolerance,  e,  a  sequence  of  step  sizes  can  be  generated, 


hi+ 2  =  q 


(hi  -  hi+ 1)  e 


I  y. 


(i+l) 

'^i+l 


(59) 


which  allows  a  numerical  ODE  solver,  such  as  those  employing  the  Dormand-Prince  method  to  minimize 
functional  error  by  adjusting  the  step  size  after  each  time  step.2  Numerical  integration  required  for  Eqs. 
19-20  and  energy  derived  from  the  circuit  model,  given  as 


were  performed  via  the  Gauss-Kronrod  quadrature  method2  at  each  time  step. 


(60) 


III.  Results 


A.  Validation  Against  Experiment 

To  validate  the  accuracy  of  the  model  described  in  this  paper,  comparisons  with  data  presented  in  Roupassov 
et.  al10  is  provided.  The  experimental  parameters  that  were  mentioned  and  used  in  the  circuit  model  are 
given  in  Table  1.  Ref.  10  uses  the  electrode  configuration  detailed  in  Fig.  8. 

Fig.  9  illustrates  an  approximation  of  the  applied  voltage  square  wave  that  was  introduced  in  the  experimen¬ 
tal  work  of  Roupassov  et.  al.10  The  slope  that  is  introduced  is  to  simulate  a  function  that  is  differentiable. 
This  is  needed  for  V^pp(t)  in  the  governing  differential  equation  as  detailed  by  Eqs.  4-5.  One  could  also 
generate  a  continuous  function  using  Fourier  decomposition  of  a  traditional  square  wave,  however  this  would 
not  account  for  the  minor  rise  and  fall  times  found  in  experimentation. 


B.  ‘Hot  Spot’  Results 

For  the  small  region  of  0.4  mm  x  0.4  mm  over  all  time  outside  of  the  anode  sheath,  the  time  variation  in 
the  number  density  can  be  established.  As  an  initial  condition  for  Eq.  28,  1015  m-3  electrons  were  assumed 
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Figure  8.  Experimental  scheme  used  by  Roupassov  et.  al10  for  the  discharge  gap.  1  high-voltage  electrode;  2 
dielectric  layer;  3  low-voltage  electrode,  4  zone  of  discharge  propagation,  5  insulating  plane. 


Table  1.  Experimental  Parameters11 


ha 

0.4  mm 

hd 

0.3  mm 

€d 

2.7 

ea 

1 

An 

30  mm2 

Ad 

30  mm2 

V 

50  kV 

AT 

65  ns 

Figure  9.  Plot  of  the  input  Voltage,  Vapp  vs.  time  used  in  the  model. 


based  on  work  by  Ref.  10.  As  displayed  in  Fig.  10,  there  is  a  large  gradient  that  occurs  during  the  rise  time 
that  peaks  around  1.13  x  1019  m-3.  It  is  also  evident  in  Fig.  10  that  the  recombination  of  electrons  is  largely 
negligible  on  the  nanosecond  time  scale.  If  a  frequency  of  1  kHz  is  used,  the  recombination  of  electrons 
with  atoms  allows  a  steady-state  electron  number  density  to  be  achieved  on  the  nanosecond  time  scale.  The 
recombination  of  electrons  becomes  a  significant  quantity  when  exploring  the  dynamics  of  a  plasma  discharge 
on  the  microsecond  time  scale. 

Using  Eq.  4  and  Fig.  5,  the  model  was  able  to  produce  a  time- varying  electron  temperature.  As  displayed 
in  Fig  11,  there  is  an  initial  spike  in  the  electron  temperature  to  29  eV  (340  000  K)  that  coincides  with 
the  peak  in  electron  number  density  at  11  ns.  Fig.  11  also  suggests  that  the  assumption  of  a  constant 
electron  temperature  is  not  an  accurate  assumption  for  nanosecond  pulsed  DBD  plasmas.  The  variability 
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Figure  10.  Plot  of  electron  number  density  vs.  time  for  the  ‘hot  spot’. 


in  the  electron  temperature  beyond  40  ns  is  due  to  the  numerical  error  tolerances  that  are  selected  when 
solving  Eqs.  4,  19  and  28.  Fig.  11  shows  that  when  reducing  the  relative  error  in  the  Dormand-Prince 
method  from  10-2  to  10-3,  the  strong  functional  variability  experiences  a  significant  reduction.  The  higher 
the  gradients  are  during  the  rise  and  fall  time,  the  finer  the  relative  error  tolerances  are  required  to  be  to 
guarantee  convergence  of  a  solution. 


- ^ - 1 - 1 - 1 - 1 - 1 - 

3 

/; 

g2'5 

: 

E  2 

Q- 

| 

c  15 
| 

" 

^  1 

_  ; 

0.5 

f  -  ,  Ou *  X~h.. .  ..4.. 

0  1  2  3  4  5  6 

Time  (s)  xi0  a 


Figure  11.  Plot  of  electron  temperature  vs.  time  for  the  ‘hot  spot’.  Left-  10  2  error.  Right-  10  3  error. 


Using  the  results  displayed  in  Figs.  10-11,  the  total  power  dissipated  to  neutral  species  as  a  function  of 
time  can  be  established.  Using  the  result  of  Eq.  4  and  the  relationship  for  the  instantaneous  power,  the 
time- varying  power  imparted  to  the  flow  can  be  given  as5 


p(t)  = 


A  V2(t) 

R/(t)  ' 


(61) 


As  shown  in  Fig.  12,  the  instantaneous  power  is  dominate  during  the  rise  time  for  the  ‘hot  spot’  region. 
Upon  integrating  the  instantaneous  power  over  time  using  Eq.  60,  this  model  produces  an  energy  value  of 
2.1  mJ  for  this  region.  When  compared  to  the  experimentally  determined  value  of  4.2  mJ  by  Roupassov  et. 
al,10  this  model  produces  an  order  of  magnitude  estimate  for  the  energy  imparted  to  neutral  species  in  this 
region. 
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Figure  12.  Plot  of  power  vs.  time  for  the  ‘hot  spot’. 


C.  ‘Hot  Spot’  Sheath  Results 

Using  the  Debye  length  approximation  provided  by  Eq.  50  and  the  solution  from  the  quasineutral  ‘hot  spot’ 
region  (ne  ~  1019  1/m3,  Te  ~  10  eV)  a  Debye  length  of  7.43  /im  is  generated.  The  electron  number  density 
experienced  in  the  anode  sheath  is  assumed  to  be  equal  to  the  values  calculated  in  the  adjacent  ‘hot  spot’ 
region.  Fig.  13  shows  the  time  variation  in  the  electron  number  density  within  the  sheath  and  also  shows  a 
peak  number  density  of  1.13  x  1019  m-3. 


Figure  13.  Plot  of  electron  number  density  vs.  time  for  the  ‘hot  spot’  sheath. 

Using  the  revised  form  of  the  electric  potential  within  the  anode  sheath  given  by  Eq.  48  and  using  Fig.  5, 
the  electric  temperature  variation  over  the  duration  of  the  pulse  can  be  established.  As  displayed  in  Fig  14, 
there  is  an  initial  spike  in  the  electron  temperature  to  (348  000  K)  that  coincides  with  the  peak  in  electron 
number  density  at  11  ns  much  like  the  hot  spot  region. 
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Figure  14.  Plot  of  electron  temperature  vs.  time  for  the  ‘hot  spot’  sheath. 


As  shown  in  Fig.  15,  the  instantaneous  power  is  dominate  during  the  rise  time  for  the  anode  sheath  region. 
Upon  integrating  the  instantaneous  power  over  time  using  Eq.  60,  this  model  produces  an  energy  value 
of  0.045  mJ  for  this  region.  This  number  is  quite  small  compared  to  the  2.1  mJ  experienced  in  the  ‘hot 
spot’  region.  However,  the  anode  sheath  does  have  a  slightly  higher  linear  energy  density  than  the  ‘hot  spot’ 
region  (6  J/m  vs.  5.25  J/m).  The  reason  that  there  is  not  significant  deviation  predicted  in  the  anode  sheath 
and  quasineutral  ‘hot  spot’  regions  is  that  explicit  charge  buildup  is  not  accounted  for  in  the  circuit  model 
within  the  sheath  region.  As  shown  by  Ref.  13,  the  electric  field  in  the  sheath  and  adjacent  quasineutral 
regions  are  approximately  equal  until  charge  buildup  is  allowed  within  the  anode  sheath  during  the  length 
of  the  pulse  or  a  series  of  pulses. 


Figure  15.  Plot  of  power  vs.  time  for  the  ‘hot  spot’  sheath. 


D.  ‘Tail’  Results 

For  the  region  0.4  mm  x  4.6  mm,  the  time  variation  in  the  number  density  can  be  established.  As  an  initial 
condition  for  Eq.  28,  1015  m-3  electrons  were  assumed,  the  same  number  of  electrons  assumed  for  the  ‘hot 
spot’  region.  When  comparing  Figs.  10  and  16,  the  tail  region  experiences  a  lower  growth  rate  in  the  number 
of  electrons  which  is  due  to  the  lower  electric  field  experienced  by  this  region.  As  described  by  Eq.  28,  a 
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lower  electric  field  produces  a  lower  number  of  ionizations  and  therefore  a  more  gradual  rise  and  lower  total 
peak  in  Ne ,  approximately  3.4  x  1017. 


Figure  16.  Plot  of  electron  number  density  vs.  time  for  the  ‘tail’. 


Using  Eq.  4  and  Fig.  5  the  model  was  able  to  produce  a  time-varying  electron  temperature.  As  displayed 
in  Fig  17,  the  traditional  assumption  of  1  eV  (11  600  K)  does  not  agree  well  with  the  results  obtained  in 
this  model  for  the  tail  region  on  the  nanosecond  time  scale.  Instead  Fig.  17  suggests  that  the  peak  electron 
temperature  is  achieved  during  the  rise  time  of  the  pulse,  220  000  K  (19  eV)  and  then  trends  downward 
during  the  plateau  portion  of  the  voltage  waveform.  Much  like  Fig.  11,  the  highest  electron  temperatures 
are  achieved  during  the  initial  high  gradient  of  the  pulse.  Fig  17.,  unlike  Fig.  11,  also  shows  a  rise  in  electron 
temperature  during  the  fall  time  of  the  pulse  as  well.  This  is  due  to  the  lower  electric  potential  and  negative 
gradient  experienced  in  the  tail  region  during  the  fall  time. 


Figure  17.  Plot  of  electron  temperature  vs.  time  for  the  ‘tail’. 
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As  shown  in  Fig.  18,  the  instantaneous  power  is  dominate  during  the  plateau  portion  of  the  applied  voltage 
pulse  for  the  ‘tail5  region.  Upon  integrating  the  instantaneous  power  over  time  using  Eq.  60,  this  model 
produces  an  energy  value  of  6.6  mJ.  When  compared  to  the  experimentally  determined  value  of  8  mJ  by 
Roupassov  et.  al,10  this  model  produces  an  absolute  error  of  ~  17.5%.  The  significant  improvement  in 
accuracy  for  the  tail  region  can  likely  be  attributed  to  the  larger  distance  from  the  cathode.  This  increase 
in  distance  improves  the  assumptions  of  quasineutrality  and  that  the  spatial  diffusion  of  charged  species  is 
negligible.  The  region  close  to  the  cathode  features  complicated  ion  and  electron  buildup  and  as  the  relative 
distance  from  the  cathode  increases,  its  impact  on  the  problem  becomes  negligible. 


Figure  18.  Plot  of  power  vs.  time  for  the  ‘tail’. 

Table  2  summarizes  the  results  obtained  using  the  circuit  model  presented  in  this  paper  and  associated 
experimental  measurements  made  by  Roupassov  et.  al10  for  both  the  ‘hot  spot’  and  ‘tail5  regions. 

Table  2.  Comparison  between  calculated  and  experimentally  measured  energy  deposition. 


Circuit  Model  [mJ] 

Experimental10  [mJ] 

Abs.  Error  [%] 

‘Hot  Spot’ 

2.1 

4.2 

50.0 

‘Tail’ 

6.6 

8.0 

17.5 

Total 

8.7 

12.2 

28.7 

IV.  Conclusion 

A  new  lumped  element  circuit  model  was  presented  that  is  valid  for  pulsed  DC  Dielectric  Barrier  Discharge 
(DBD)  plasmas.  The  model  approximates  the  total  energy  dissipated  into  neutral  species  using  a  lumped 
element  circuit  while  containing  relevant  plasma  physics  in  the  form  of  a  variable  electron  temperature 
and  number  density.  An  approximate  expression  was  formulated  using  the  conductivity  of  the  discharge 
to  calculate  the  resistance  value  for  the  air  gap.  Asymmetric  wall  effects  were  also  approximated  in  the 
model  by  including  the  effect  of  the  anode  sheath.  Results  of  the  model  were  verified  against  a  pulsed  DC 
experiment  conducted  by  Roupassov  et.  al10  and  order  of  magnitude  agreement  was  obtained  for  the  energy 
imparted  into  the  plasma  in  both  the  homogeneous  ‘hot  spot’  region  and  ‘tail5  region. 

Appendix 

The  momentum  cross  sections  used  for  the  numerical  approximations  in  the  model  are  shown  in  Table  3. 
These  values  were  obtained  by  A.V.  Phelps.8 
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Table  3.  Momentum  Cross  Sections8 


Kinetic  Energy  (eV) 

N2  (IE-16  cm2) 

O2  (IE-16  cm2) 

1 

10 

7.2 

2.1 

27 

6.6 

3 

21.7 

5.7 

4 

12.6 

5.5 

5 

10.9 

5.6 

10 

10.4 

5 

15 

11 

8.8 

20 

10.2 

8.6 
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